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Abstract 

Nonlinear stages of the recently uncovered instability due to insoluble surfactant 
at the interface between two fluids are investigated for the case of a creeping plane 
Couette flow with one of the fluids a thin film and the other one a much thicker layer. 
Numerical simulation of strongly nonlinear longwave evolution equations which couple 
the film thickness and the surfactant concentration reveals that in contrast to all similar 
instabilities of surfactant-free flows, no amount of the interfacial shear rate can lead to a 
small-amplitude saturation of the instability. Thus, the flow is stable when the shear is 
zero, but with non-zero shear rates, no matter how small or large (while remaining below 
an upper limit set by the assumption of creeping flow), it will reach large deviations 
from the base values- of the order of the latter or larger, ft is conjectured that the 
time this evolution takes grows to infinity as the interfacial shear approaches zero. It is 
verified that the absence of small-amplitude saturation is not a singularity of the zero 
surface diffusivity of the interfacial surfactant. 
AMS Subject Classification: 76E17,76E30 
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1. Introduction 



Recently, a new instability of a channel flow of two fluid layers was uncovered in Frenkel 
& Halpern Halpern & Frenkel Q. This linear instability was further investigated in 
Blyth & Pozrikidis B, Blyth & Pozrikidis 0], Wei 0, Frenkel & Halpern 0, Wei 
(Blyth & Pozrikidis ja| and Blyth & Pozrikidis included some nonlinear investigations). 
It is driven by the interaction of capillary and Marangoni forces with the interfacial shear of 



the base flow, and does not depend on gravity and inertia for its realization. The instability 
disappears if the shear is zero. On the other hand, as was noted in Frenkel & Halpern 
[]J, a nonzero shear is known to enable the nonlinear saturation of different instabilities, 
in which capillarity plays a stabilizing role and the destabilizing factor can be molecular 

Van der Waals forces (Babchin et al or gravity (Babchin et al M, Halpern & Frenkel 

or the interfacial jump in viscosity (Hooper & Grimshaw [ill ]! Shlang et al 0])> 
or capillary forces due to the azimuthal curvature of the interface in core-annular Hows 
(Frenkel et al [l3j], Papageorgiou et al [lI. Halpern & Grotberg |l5|) (see also reviews 
Joseph & Renardy [16J, Joseph et al [lTj]). When the shear rate is sufficiently large, the 
instabilities saturate with the final amplitude of the interface undulations being small, so 
that the whole evolution starting from small disturbances can be described by a weakly- 
nonlinear equation for the film thickness. (The mechanism of this nonlinear saturation was 
proposed in Babchin et al [j| and further clarified in Frenkel et al In other cases, 

the saturated amplitudes for certain regimes are large and a strongly nonlinear equation is 
appropriate, as for example for the flow down a fiber (Frenkel 0]; Frenkel [l^, Kerchman 
& Frenkel j^, Kalliadasis & Chang |2j|, Chang & Dcmckhin 22], Kliakhandler et al 
2^) or for the core-annular flow (Frenkel & Kerchman j^], Kerchman j^, Joseph et al 
171]). However, for the same systems, there are other parametric regimes for which small- 
amplitude saturation still occurs, so that the weakly-nonlinear description is good. Since 
for the instability considered below the interfacial shear of velocity plays two conflicting 
roles, working both for the linear instability and the nonlinear mechanism of saturation 
(and also because the dynamics involves a surfactant evolution equation in addition to the 
film thickness equation, in contrast to the single-equation description of evolution in the 
surfactant-free instabilities), there is a question whether there are any parametric domains 
for which the small-amplitude saturation of the interfacial-surfactant instability happens 
(Frenkel & Halpern Q). In this paper, we investigate the nonlinear stages of the instability 
with respect to this question (and, to anticipate, arrive at the conclusion that no small- 
amplitude saturation occurs, so that the weakly-nonlinear description inevitably breaks down 
after some time). 

The rest of the paper is organized as follows. In the next section, after giving the exact 
formulation of the Navier-Stokes problem, we use the lubrication approximation to obtain 
a simplified system of coupled equations for the film thickness and the surfactant concen- 
tration. The linear stability governed by these equations is studied in section 3 and the 
nonlinear evolution is numerically simulated in section 4. The last section contains the 
discussion of results and conclusions. 



2. Problem formulation 



2.1 Exact dynamic equations 

The formulation of the exact Navier-Stokes problem (as given earlier in Frenkel & Halpern 
0, Halpern & Frenkel 0]) is as follows. Consider two immiscible fluid layers between two 
infinite parallel plates, as in Fig. ^ Let the base flow be driven by the combined action of 
an in-plane steady motion of one of the plates and a constant pressure gradient parallel to 
the plate velocity and directed in either the same or the opposite sense. It is well known 
that the base "Couette-Poiseuille" velocity profiles are steady and vary (quadratically) in 
the spanwise direction only, and the base interface between the fluids is flat. For simplicity, 
let the densities of the two fluids be equal. Then gravity does not affect stability of the base 
flow, and is disregarded below. Let y* be the spanwise, "vertical" , coordinate (the symbol * 
indicates a dimensional quantity). Let the interface be at y* — di where di is the thickness 




Figure 1: Definition sketch for a viscous two-layer Couette flow between a fixed plate at 
y* = and a moving plate at y* = d\ + d 2 . The disturbed interface y* = h*(x* , V) is shown 
as a sinusoidal curve. T*(x*,t*) is the concentration of the insoluble surfactant monolayer 
shown as discs with tails. The base velocity profiles are shown by bold arrows. 



of the thinner layer, and the y*-axis directed from the thinner layer ("film") to the thicker 
one; we will call this the "upward" direction (clearly, since there is no gravity, the notions of 
"up" and "down" are a matter of convention). Thus, d\ < d 2 holds, where d 2 is the thickness 
of the upper fluid. The direction of the "horizontal" £*-axis is chosen so that velocity at the 
interface is positive, say V\ (whereas the velocity at the lower plate is zero). 

The Squire-type theorem for the case with surfactants, proved in Halpern & Frenkel Q, 
allows us to confine the stability considerations to two-dimensional perturbed flows (in the 
£*j/*-plane). The equation of the interface is y* = h*(x* , t*), where h* is the film thickness. 
The Navier-Stokes and incompressibility equations governing the fluid motion in the two 
layers are (with j — I for the lower layer and j = 2 for the upper one) 



<9v* 



v* • V*v*) = -V*p* 



V* • v* = 



where V* 



_d 

^ dx* ' dy 



t), p is the density (of both fluids), v* = (u*,v*) is the fluid velocity 



with the horizontal component u* 



and vertical component v* , and p* is the pressure. 



We use the "no-slip , no-penetration" boundary conditions (requiring zero relative velocities) 
at the plates: u\ = 0, v\ = at y* = 0; and u\ = U 2 , v\ = at y* = d\ + d 2 , where U 2 is 
the velocity of the upper plate (which can be positive or negative for the general Couette- 
Poiseuille flow). The interfacial boundary conditions are as follows. The velocity must be 
continuous at the interface: [v*]f = 0, where [A\\ = A 2 — A\ denotes the jump in A across 
the interface, i.e. at y* — h*(x*,t*). The interfacial balances of the tangential and normal 
stresses taking into account the jump in the tangential stress (due to the variability of surface 



tension) as well as the capillary jump in the normal stress, at y* — h*(x* , t*), are 



' -[(1 - K 2 ,Uu* + „*.) + 2hf^{vl, - <,)]? 



l + /i*? LV I s * ; » xm (I + /1*?) 1 / 2 

[(i + h%y 2/x(^<, - /£.( u ;, + <.) + «;.)]? = (i+i'*V)V2 g * 

where a* is the surface tension (the subscripts x* , y* , f * denote the corresponding differen- 
tiation) . The kinematic boundary condition (the conservation of mass condition) is 



h'(x* ,t") 

dh* d 
dt* dx* 



u*dy* = 0. (1) 



The surface concentration of the insoluble surfactant on the interface, T*(x* , t*), obeys the 
following equation (see a simple derivation in Halpern & Frenkel 0): 

^ r *> » (a rv,_ „» (1£) W 



<9i* 9x* 9a;* \ H dx 



where = yl-j- /i* 2 , and D s is the surface molecular diffusivity of surfactant; D s is usually 
negligible, and is discarded below. We assume that the surfactant concentration T*(x*,t*) 
is sufficiently small, remaining always far below its saturation value, so that the linear 
dependence of the surface tension on the surfactant concentration is a good approximation 
(see Edwards et al p. 171): a* = (To — E (T* — To), where To is the base surfactant 

concentration, ao is the base surface tension and E is a constant. 

We use the following nondimcnsionalization: 

(x*,y*) t* (u*,v*) p* T* a* 

{x, y) = — -, — , t = j— , {u, v) = , p = — 7j- , r = — , a = — . (3) 



In the base flow, the interface is flat, h = 1, and the surfactant concentration is uniform, 
T = 1, where the overbar indicates a base-state quantity. Since the base flow has a constant 
pressure gradient, say 2q, the two layers have equal pressures, p 1 = p 2 = 2qx, and the 
velocity profiles are 

My) =ry + qy 2 , vi=0, forQ<y<l, (4) 

ff) \ f Q 

u 2 (y) = (r + q) + —y+—y 2 , v 2 =0, forl<y<n+l, (5) 

m m m 

where n = d 2 jd\ is the aspect ratio, n > 1, and m — ^2/^1- The pair of constants r and q 
is used in place of the physical control parameters, the relative velocity of the plates and the 
pressure gradient, to characterize the base flow. The base interfacial velocity is w = r + q, 
so r + q > 0. The parameter s = r + 2q is the interfacial shear of velocity (on the lower-fluid 
side). 

2.2 Longwave (lubrication) approximation 

We will assume that the aspect ratio is large, n >> 1 (actually we work in the zeroth order of 
the small parameter l/n, so the parameter n disappears from the equations). For sufficiently 
slow flow, the inertia terms in the Navier-Stokes equations are negligible. In the well-known 
lubrication approximation, which implies that the streamwise characteristic lengthscale is 



much larger than the film thickness (see, e.g., the review papers Frenkel & Indircshkumar 
|27j and Oron et al |28j). the simplified dynamic equations are, in the dimensionless form, 



dp £Pu 

dx 8y 2 ' 1 ' 

« 

dv du 
dy dx' 

where u = u\ — Ui, v = Vi — Vi, and p = pi ~Pi- 

The no-slip no-penetration boundary conditions at the lower plate prescribe that 

u = = v at y = 0. 

As usual for this class of systems with infinitel y la rge aspect ratio (e.g., Babchin et al 
Babchin et al 9], Hooper [2^|, Frenkel et al Jlj||). the disturbance quantities of the thick 
layer are negligible in the interfacial stress boundary conditions (at least if the viscosity 
ratio is not too different from one, so that at y = h(x, t) the difference of base velocities 
is not much larger than u, and hence u and the corresponding velocity disturbance in the 
thick layer are of the same order of magnitude), so the simplified boundary conditions at 
y = h{x, t) for the system 110181) are 

p = -[l-M{T-l)]h xx 

(the curvature term is retained to allow for the possibility of large surface values of tension) 
and 

u y = -MT X 

where the Marangoni number M is defined as M = ETq/o-q, and we have neglected h x (as 
compared to 1) as well as some other terms, in accordance with the longwave (small-slope) 
approximation. The solution for u in terms of h(x,t) and T(x,t) is 

u = {h[(l + M- MT)h xx ] x - MT x }y - + M - MT)h xx ] x y 2 . 

We substitute the total velocity u\ — ry+qy 2 +u into (i) the dimensionless mass conservation 

h(x,t) 

equation ifTjl. ht + [ J U\(x, y, t) dy] x — 0, and (ii) the simplified surfactant conservation 
o 

equation 

T t + [T{rh + qh 2 + u)] x = 

(with u being evaluated at the interface y = h(x,t)), where we have set H = 1 (since 
h 2 << 1) and the right-hand side of equation J3J) to zero since the molecular diffusivity D s 
is usually very small. As a result, we obtain a system of two coupled evolution equations, 
for h and T: 

h t + (Vy +gy-Mr x y + ^[(1 + M - MT)h xx } x hA =0, (9) 



r 



•h + qh 2 - MT x h + hi 2 [(l + M - MT)h xx l 



= 0. (10) 



The terms containing r and q are clearly due to the nonzero base velocity. 



To simplify, we make some additional assumptions: We assume q — (then the velocity 
is piecewise-linear and r is its slope in the film; also r is equal to the intcrfacial velocity) 
and M << 1, and also MT « 1. Rescaling x = fix and t = r(3t where (3 = (3M/2) 1 / 2 , 
introducing the constant C — /3M/2r = (3M 3 /(8r 2 )) 1//2 , and dropping the tildes from the 
new variables, the equations take the form 

h 2 

ht + [ y + C(-T x h 2 + h xxx h 3 )] x = 0, (11) 



r t + 



r ( h + c(-2T x h + h 2 h xxx ) 



0. (12) 



In each of these equations, the second term is due to the interfacial shear, the third one is 
due to Marangoni forces, and the last one describes the effect of capillarity. 

In the rescaled variables, the longwave requirement reads 

\Ph x \ « 1 

which, assuming 

Mi/2 <<cl 

can be satisfied even in some cases with h x > 1 (note that the 
unit of length for x is much greater than that for h) . 

Note that because of the conservative form of the equations l|ll|) and (|12fl . the integral 
quantities J h(x,t)dx and J T(x,t)dx are time-independent if they are taken over domains 
with periodic boundary conditions. 

We note that it is not difficult to obtain weakly-nonlinear evolution equations for the dis- 
turbances r\ and g defined by 

h=l + T), T = l+g. 
One assumes the disturbances to be finite but small, 

i] « 1, g « 1, 

and retains only the leading-order linear and nonlinear terms in r\ and g in equations (|11|1 
and (|12|) . Changing to the coordinate z — x — t (which eliminates the term r\ x from equation 
<|11|) and the term g x from equation 1121 ) . we obtain 

Vt + (y - °9* + Cr lzz^j = 0, (13) 

9t+(v9 + V- 2Cg z + ^-ijzzz^j = 0. (14) 

(In this derivation, for example the term proportional to g z r] z coming from the third term in 
Eq. (|11|) is smaller by a factor of r\ << 1 than the term proportional to g zz arising from the 
same term in Eq. i|ll|) . Therefore, the smaller term has been dropped, while, on the other 
hand, the term with (r] 2 ) z has been retained in Eq. I|13|l as the leading order nonlinear term 
since there is no other term which would be clearly larger than it: Indeed, depending on the 
given conditions, the parameter C can take values ranging from very small to very large, and 
the same is true of the characteristic length scale measured in the units we have chosen. The 
only restriction on C is imposed by the requirement of smallness for the modified Reynolds 
number, Red\/K << 1, where Re is the Reynolds number based on the interfacial velocity 
and the thickness di, and A/d± is the dimensionless streamwise length scale of solutions. 



This gives the constraint C » Af 3 / 2 pood 2 / '(/if A). The nonlinear term in Eq. H14() is clearly 
smaller than the linear term r\ z in view of g « 1; it could be omitted and has only been 
retained in order to show the leading nonlinear term.) 

These weakly-nonlinear evolution equations can be rescaled to a "canonical form" having 
just one parametric coefficient: 

r\t - g zz + Vzzzz = 0, (15) 

3 

gt ~ g zz + -^t)zzzz + Ariz = (16) 

(where the parameter A is inversely proportional to the parameter C) . [We notice that these 
weakly-nonlinear equations have certain similarities with the model system constructed for 
specific purposes in Malomed et al [3(| : Both systems are first order in time and fourth order 
in space and have similar nonlinear terms, although couplings are different. For a review of 
other systems of coupled evolution equations, see, e.g., Glasman et al However, the 

above weakly-nonlinear equations turn out to be good for just a short transient stage: as we 
will see below, the small-amplitude saturation does not take place and disturbances quickly 
grow beyond the scope of the weakly-nonlinear equations. To repeat, this absence of the 
small-amplitude saturation in the interfacial-surfactant instability, contrasting with the pres- 
ence of such small-amplitude saturation in all previously known surfactant-free instabilities 
of two-fluid flows, is the main point of the present communication. 



3. Normal modes of infinitesimal disturbances 



To examine behavior of infinitesimal disturbances to the steady uniform flow with T = 1 
and h = 1, we substitute h = 1 + r\ and T = 1 + g into Eqs. I|ll|l - (|12|l , and neglect the 
nonlinear terms, obtaining the linear system 

rjt + rix-Cg xx + Cr] xxxx = 0, (17) 
3 

gt + gx + Vx- 2Cg xx + -Crj xxxx = 0. (18) 

It is known (Frenkel & Halpern Q], Halpern & Frenkel Q) that the normal modes (77, (?) oc 
e iax e ~ft e iujt w fa eie a j s wavenumber, 7(a) is the growth rate, and uj(a) is the time 
frequency, are unstable-if and only if C < 00 (i.e., the shear velocity r ^ 0)-for "long" waves, 
< a < ao, where ao is a marginal wavenumber. The longwave asymptotic dependence of 
the growth rate is 7 oc a 3 / 2 as a — > 0. The normal modes decay if a > ao, with 7 ~ — jet 2 
as a — > 00. (We note that the decay is weaker than a 4 despite the presence of the fourth- 
derivative terms in the linearized equations l|17fl and (|18fl ). Typical dispersion curves, the 
growth rate versus wavenumber, are shown in Fig. [2J They have the maximum value j m ax 
at a certain wavenumber a max depending on the value of the equation parameter C . These 
dependences are shown in Fig. |21 

Multiplying Eq. i|17[l by 7/ and (|18|) and integrating over the interval of spatial periodicity, and 
taking into account that r\ and g are out of phase by nearly 7r in unstable modes, shows that 
the fourth derivative (capillary) term is stabilizing in Eq. (|17fl but destabilizing in Eq. i|18|) . 
Similarly, the Marangoni term, the one with g xx , is stabilizing for g but destabilizing with 
regard to r/. 
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Figure 2: Dependence of the growth rate 7 on the wavenumber a of a normal mode for three 
values of the parameter C indicated in the legend. 
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Figure 3: (a) Growth rate and (b) wavenumber of the fastest-growing normal mode versus 
C. 



4. Nonlinear numerical simulations 



We solved the system of equations and H12|l numerically on < x < 2ttQ with periodic 
boundary conditions using the method of lines. This involved the approximation of the 
spatial derivatives with finite differences and an implicit Gear's method for the time-stepping. 
The number of nodes in the spatial discretization ranged from hundreds to tens of thousands. 
Also, some results were checked using standard pseudospectral Fourier method (see e.g. 
Kerchman & Frenkel 20] ) in conjunction with an implicit Gear's method for the time- 
stepping. The number of modes ranged from 256 to 1024 for different cases. We have used 
dealiasing as described e.g. in Kerchman & Frenkel [2(J (similar to Canuto et al |32^h 
filtering out slightly more than three-fifths of the modes with higher wavenumbers from the 
unknown functions and their derivatives, computed in the Fourier space, before computing 
the nonlinear terms of the equations in the physical space. In all cases, we took care that the 
length of the interval determined by the parameter Q was sufficiently large to accommodate 
five or more elementary pulses, as it is known, e.g. for the Kuramoto-Sivashinsky equation, 
that the results can be sensitive to the type of boundary conditions and the length of the 
interval when the latter is not sufficiently large (Papageorgiou and Smyrlis 33] , Frenkel & 
Indireshkumar |2]]], Wittenber & Holmes [13] )• [ m particular, the preliminary nonlinear 
saturation results found-although for very different parametric regimes-in Blyth & Pozrikidis 
Q, may not hold for sufficiently large domains. We have checked that such was indeed the 
case for our semi-infinite flow configuration: The small-amplitude saturation occurring for a 
narrow band of values of Q near criticality, such that the periodicity interval accommodated 
just one unstable normal mode, disappeared for larger values of Q, since the corresponding 
larger intervals can accommodate additional longwave modes which are unstable. Of course, 
to any one-pulse solution for a small Q = a there is a corresponding solution for an n times 
larger value, Q = na, consisting of n identical a-pulses. However, as a rule, this n-pulse 
solution is unstable to longwave disturbances and therefore is not observed as an outcome of 
evolution. Conversely, the results are essentially insensitive to further increase of Q chosen 
in this way] 

For the initial conditions we use small-amplitude white-noise deviations from the base uni- 
form profiles. The deviation values at the spatial nodes are chosen from the uniform random 
number distribution on the interval [— 10 -2 , I0~ 2 ]. Typical initial profiles are shown in Fig.^J 

In the subsequent evolution, the shortwave normal modes die out, but the longwave ones 
grow, and the spatial profiles evolve to the pulse trains, like the travelling waves shown in 
Fig. Efa). A typical evolution of the instantaneous global maxima and minima of the spatial 
profiles of the surfactant concentration T and the film thickness h is shown in Fig. [SJbc): 
There, the deviations of T and h saturate but at levels which are not small as compared 
to non-disturbed values. In Fig. |HJ the evolution for different parameter values is shown 
to emphasize the point that the deviations in all cases grow to levels of order one. To 
this end, each time series is followed up to the instant at which the minimum surfactant 
concentration reaches the representative level of 0.5. We used five different values of the 
equation parameter C, from small to order one to large: 10 -2 , 10 _1 , 1, 10, and 10 2 , with 
sufficiently large values of Q (as discussed above; also, in all cases, we have checked that 
refining the spatial and temporal resolutions does not change the results.). The figure clearly 
shows that in all cases the deviations of the surfactant concentration and the film thickness 
grow and become non-small in comparison with their base values. As a result, the weakly 
nonlinear equations (|I5|I and (|lfci|) . based on the assumption that the deviations are small, 
would cease to be good. (Also, clearly, if the film thickness becomes sufficiently small as for 
C = 100 (see Fig. E^b)), the molecular van der Waals forces should be taken into account.) 
The greater the value of C, the larger time (of the governing equations) the evolution takes. 
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Figure 4: Typical initial profile of surfactant concentration or film thickness: a white noise 
distribution, with deviation amplitude of 10~ 2 for 512 modes. 

As was mentioned above, the initial configuration with many maxima and minima (see 
Fig. H} quickly evolves, by coalescences, to just a few pulses (see Fig. Et)- Another way to 
describe this is to note that the normal modes with large wavenumbers are damped and die 
out. The normal mode of the fastest growth has a tendency to dominate, at least initially. 
We can define an effective wavenumber cx e tf of an instantaneous profile to be 2ir over the 
"length per one pulse" . The latter is the interval length 2ttQ divided by the number of pulses 
N. Therefore, a e // = N/Q. The dependence of a e ff on C is shown in Fig. Ufa) together 
with that of the wavenumber a ma x of the fastest growing normal mode. These are seen to 
be close throughout the entire range of C shown in Fig. Ufa). 

We also define an effective growth rate j e f / to be that of a (hypothetical) normal mode that 
would have the same minimum surfactant concentration as the actual nonlinear profile at 
two instances of time, t = and t = t\, where t\ is determined by T m i n (ti) — 0.5. This 
gives j e ff = log(0.5/10 -2 )/ti, that is j e ff = log(50)/fi. Fig. [7{b) shows j e fj (inversely 
proportional to the evolution time ii) as a function of C . For comparison, the growth rate 
Imax of the fastest growing normal mode is included in the figure. Jeff is always smaller 
than j m ax because the initial exponential growth of small disturbances slows down as the 
magnitude of the disturbances grows. The figure also reflects on the fact that the evolution 
time grows with C. We conjecture that t\ — > oo as C — > oo- that is, (assuming that M is 
fixed), as the shear rate approaches zero. (It is easy to see that this would mean that the 
corresponding physical, dimensional, time of the evolution diverges to infinity as well.) 

The above results have been obtained discarding surfactant diffusion. To check the continuity 
in the surfactant diffusivity at its zero value, we should add into the left-hand side of equation 
H12(l the term — DT XX) where D is a constant proportional to D s , the surface molecular 
diffusivity of surfactant. We have simulated the modified equations with the values of D as 
large as 0.1. The results turned out to be qualitatively the same as with D = (see Fig.|SJ). 

5. Discussion and conclusions 

We have investigated the nonlinear stages of the insoluble-surfactant instability in a sheared 
viscous film in a semi-infinite plane-Couette flow. Numerical simulations of strongly non- 
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Figure 5: A typical evolution of the surfactant concentration and film thickness, (a) The 
saturated travelling-wave profiles of the surfactant concentration and the film thickness, 
at t = 5000. (b) The time dependence of the minimum and maximum (over x) surfactant 
concentration, (c) The time dependence of the minimum and maximum film thickness. Here 
C = 10~ 2 (and Q = 5.2). 




Figure 6: Evolution of (a) surfactant concentration and (b) film thickness: The dependence 
of the minimum and maximum (over x) values T max , T m i n and h maxi h m i n on time t, for 
the five different values of C shown in the legend. Each evolution is shown up to the time 
at which the minimum surfactant concentration reaches the value of 0.5. (In the order of 
increasing C, the values of Q used in these simulations are 5.2, 8, 18, 50, 250.) 
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Figure 7: C-dependence of (a) the effective wavenumber and (b) the effective growth rate 
(defined in the text). For comparison the corresponding linear-theory quantities ^max and 
oimax are included. 
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Figure 8: Same as Fig. 6 but for a non-zero diffusivity D value. Here D 
(Q = 18). 



0.1, C = 1 



linear longwave evolution equations on extended spatial intervals show that the earlier un- 
covered remarkable property of the linear instability, the necessity of a non-zero interfacial 
shear rate for the system to be unstable to infinitesimal disturbances, has a "nonlinear con- 
tinuation" : no matter how large or small the interfacial shear of the base flow, there is no 
small-amplitude saturation of the instability. (So the weakly non-linear equations cannot 
describe the large-time evolution.) 

The impossibility of a small-amplitude saturation is corroborated by the following argument. 
Suppose there is such a saturated state of Eqs. (|13|l and 1(14(1 with N << 1 and G « 1, 
where N and G are the saturated amplitudes of deviations of the film thickness and surfactant 
concentration, respectively, having the characteristic length scale L (such that r\ z ~ N/L 
and g z ~ G/L) and time scale T (such that rjt ~ N/T and g t ~ G/T). Consider the 
dominant balance of terms in Eqs. ((13(1 and (|14(l . The only nonlinear term, rjf] z , must be 
one of the dominant terms in Eq. (13) ■ But at least one of the terms Cg zz and Cr\ zzzz 
must be among the dominant terms there, since the dominant equation cannot be just 
f]t + Viz = 0: this equation is well-known to lead to infinite slopes in finite time. Thus, 
N 2 /L ~ max(CG / L 2 , CN/L 4 ). But the same terms, Cg zz and Ci] zzzz , cannot be among the 
dominant ones in Eq. 1)14(1. since the term r\ z there is much larger: N/L » N 2 /L. Thus the 
dominant balance in Eq. 1)14(1 would have to be just gt+r]z = 0, leading to r\ z ~ gt- Combining 
this relation with r\ t ^ VVz ( see Eq. (fT3|) 1. we can write rjt < r\gt ~ NG/T « N/T ~ r\ t . 
Thus, we have arrived at rjt « ??t, a contradiction. Therefore, there can be no saturated 
solutions of the weakly- nonlinear equations ((13(1 - 1(14(1 with both n and g being small; in 
other words, there can be no small-amplitude saturation of the instability in question. (This 
argument assumes a single time scale, which is generically the case. Only for a limited 
domain of near-critical states multiple time scales occur due to the presence of the small 
super-criticality parameter. As was mentioned above, under such non generic conditions the 
small-amplitude saturation does occur. This should be possible be describable by Landau 
type theory, which we plan to investigate in the future.) 



The same considerations hold for the weakly nonlinear equations obtainable directly from 



Eqs. PT71|) . They imply that, even when the Marangoni number M is not small, the small- 
amplitude saturation of the interfacial-surfactant instability is impossible as well. 

Our numerics show that when one includes into consideration (small) non-zero diffusion of 
the surfactant, the instability still does not saturate with small amplitudes, and the distur- 
bances of the film thickness or/and the surfactant concentration reach values of order one. 
Evidently, the increase of the base shear does not lead to the small-amplitude saturation-as 
it does in all previously known similar systems without surfactants (which are described by a 
single evolution equation, the one for the film thickness). Thus, for the interfacial surfactant 
instability under consideration here, the increase of the saturating effect with the shear rate 
cannot overcome the simultaneous increase in the strength of the instability. 

It remains to be seen if there is a threshold in the strength of the surfactant diffusivity above 
which a small-amplitude saturation of the instability will set in. We hope to fully investigate 
the effect of diffusion (including the questions of saturation, with small as well as with large 
amplitudes) in the future. Even without diffusion, the large-amplitude saturation is possible; 
an example is shown in Fig[5] The systematic investigation of the large-amplitude evolution, 
which is beyond the scope of this paper, will be a subject of the future study as well. 
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